Mon. Not. R. Astron. Soc. 000, 000,000 (0000) 



Printed 1 February 2008 



(MN I^TeX style file vl.4) 



Modelling the evolution of correlation functions in 
gravitational clustering 

Dipak Munshi 1 , T. Padmanabhan 2 

1 Queen Mary and Westfield College, London El 4NS, United Kingdom 

2 Inter-University Centre for Astronomy and Astrophysics, Post Bag 4, Ganeshkhind, Pune 411 007, India 
Email: D.Munshi@qmw.ac.uk, nabhan(3iucaa.ernet.in 



1 February 2008 



Key words: 

analytical 



On 

m 

> 
O 

o 

Oh 

O ' 1 INTRODUCTION 

a: 

There is growing evidence that the large scale structure in 
r* ■ the universe formed through gravitational amplification of 
. £h ' sma U inhomogeneities. Semianalytic modelling of gravita- 
, tional clustering of collisionless, non relativistic, dark mat- 
5_l ■ ter particles will be of significant utility in understanding 
CS ' the formation of large scale structures. Such a modelling is 
straightforward when the density contrasts are small and 
perturbation theory, based on a suitably chosen small pa- 
rameter, is applicable (see e.g., Fry 1984, Moutarde et al 
1991, Buchert 1992, Bernardeau 1992 ). In the other ex- 
treme, highly nonlinear regimes can be handled if one is 
prepared to make some extra assumptions like stable clus- 
tering ( Peebles 1980 ) or those which underly the Press- 
Schecter formalism, Peaks formalism etc. ( Bardeen et al. 
1986 ). The intermediate regime is considerably more diffi- 
cult but some progress has been made recently even in this 
case ( Padmanabhan, 1996; also see Hamilton et al. 1991, 
Nityananda and Padmanahban 1994), using the scale invari- 
ant spherical infall models. These papers give an expression 
for the nonlinear mean correlation function in terms of the 
linear mean correlation function both in the intermediate 
and nonlinear regimes. To do so, Padmanabhan (1996) has 



ABSTRACT 

Padmanabhan (1996) has suggested a model to relate the nonlinear two - point corre- 
lation function to the linear two - point correlation function. In this paper, we extend 
this model in two directions: (1) By averaging over the initial Gaussian distribution 
of density contrasts, we estimate the spectral dependence of the scaling between non- 
linear and linear correlation functions. (2) By using a physically motivated ansatz, we 
generalise the model to N-point correlation functions and relate the nonlinear, volume 
averaged, N-point correlation function £,n(x, a) with linearly extrapolated volume av- 
eraged 2-point correlation function a) evaluated at a different scale. We compare 
the point of transition between different regimes obtained from our model with nu- 
merical simulations and show that the spectral dependece of the scaling relations seen 
in the simulations can be easily understood. Comparison of the calculated form of £jv 
with the simulations show reasonable agreement. We discuss several implications of 
the results. 
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concentrated on a typical spherical region and has ignored 
the effects arising out of averaging over peaks of different 
sizes. Also no attempt was made to model higher order cor- 
relation functions. In this paper, we generalise these ideas 
in two directions: 

(a) We consider peaks of different heights and average 
over them along the lines suggested in Padmanabhan et al. 
(1996). Such an averaging will introduce spectrum depen- 
dence in the relation connecting nonlinear and linear corre- 
lation functions which have been noted in numerical simula- 
tions ( Padmanabhan et al. 1996, Jain et al. 1995, Peacock 
& Dodds 1996). We will show that this dependence can be 
understood from our model. 

(b) We shall postulate an ansatz for the higher order 
correlation functions, generalising the result for two point 
correlation function. Using this ansatz, we shall calculate 
the Sjv parameters for both the intermediate and nonlinear 
regimes and compare them with the simulations. We will 
see that there is reasonable - though not excellent - agree- 
ment between the model and simulations suggesting that 
the ansatz we have proposed is along the right direction. 
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2 THE MODEL AND THE ANSATZ 

The basic idea behind the model used in Padmanabhan ( 
1996 ) can be described as follows: Consider the evolution of 
density perturbations starting from an initial configuration 
which we take to be a realisation of a Gaussian random field 
with variance a. A region with initial density contrast Si will 
expand to a maximum radius x ta = Xi/$i and will finally 
collapse to an object of radius Xf which will contribute to 
the two-point correlation function an amount proportional 
to (xi/xf) 3 . The initial density contrast within a randomly 
placed sphere of radius Xi will be vo{xi) with a probability 
proportional to exp(— v 2 j2). On the other hand, the ini- 
tial density contrast within a sphere of radius Xi, centered 
around a peak in the density field will be proportional to the 
two-point correlation function and will be v 2 ^{xi) with a 
probability proportional to exp(— v 2 /2). [We have obtained 
the quadratic scaling in u based on the assumption that £ 
scales in the same way as mean square fluctuations in the 
mass, which - in turn - will scale as the mean square of the 
gaussian density field. In general, one expects the scaling to 
be u a with a « 2 — 1. The results are easily generalisable 
to any value of a. We will stick to a = 2 since it gives rea- 
sonable agreement with simulations and is based on simple 
considerations ] . It follows that the contribution from a typ- 
ical region will scale as £ ni oc t 3 ^ 2 while that from higher 
peaks will scale as £„j oc £f . In the intermediate phase, most 
dominant contribution arises from high peaks and we find 
the scaling to be £ n ; oc £f. The non- linear virialized regime 
is dominated by contribution from several typical initial re- 
gions and has the scaling £ ni oc £^ 2 . This was essentially the 
feature pointed out in Padmanabhan (1996) though in that 
work it was assumed that v — 1. To take into account the 
statistical fluctuations of the initial Gaussian field we can 
average over different v with a Gaussian probability dis- 
tribution. [Strictly speaking, there will be deviations from 
pure gaussian distribution because our averaging requires a 
mapping from lagrangian to eulerian coordinates; we shall 
ignore this because it is a higher order effect]. We shall do 
this calculation in the next section. 

To generalise the above ideas for higher order corre- 
lation functions is more nontrivial. In general n-point cor- 
relation functions will depend on shapes but volume aver- 
aging will remove this shape dependence. The "Sn param- 
eters" are then defined as dimensionless ratios of £jv(:r,a) 
and &(x, a)' iv ~ 1 ' ) . Such volume-averaged N-point functions 
( which can be directly related to counts-in-cells ) and the 
Sn parameters have been studied extensively in literature ( 
White 1979, Balian & Schaeffer 1984, Bouchet et al. 1991, 
Bouchet & Hernquist 1992) . The Sn parameters show fairly 
simple pattern of behaviour both in the perturbative and 
nonlinear regimes. It can be shown that all Sn's can be 
evaluated from spherical collapse model in the limit £2 — * 0. 
In this limit, they are constant and depend only on initial 
spectral index when smoothing is taken into account. They 
are also expected to be constants in the nonlinear regime. 
These results indicate that the hierarchical pattern, which is 
generally assumed to describe nonlinear £jv functions, could 
have a larger range of validity. We shall exploit this possi- 
bility to estimate the Sn in the intermediate and nonlinear 
regimes along the following lines: 

The evolution of N-point correlation functions is de- 



scribed by momentum moments of BBGKY hierarchical 
equations, which can be expressed in the form 

Here a varies from 1 to N, i varies over the Cartesian compo- 
nents and Qn is the full N-point correlation function given 
by 

Q N {1,2,...N) = 1 + cT 2 (l,2) + .... 

+ 6(1,2,3) + ....+ £jv(1,2,3,...AT) (2) 

and £jv denotes the reduced part of N-point correlation func- 
tion. For 2-point correlation function, the resulting equation 
can be simplified to ( Peebles 1980 ) 

which describes the conservation of pairs. In the integral 
form, the same result can be expressed as 

x 3 (l + &(x)) = l 3 (4) 

where I = {x 3 ) 1 ^ 3 is the average initial scale from which 
collapsed structures of size x have formed. 

Our aim is to generalise the above result for higher or- 
der correlation functions, but it is obvious that one can not 
get such a simple relation for higher order correlation func- 
tions in which N — 1 different length scales are present. To 
make progress, one needs to assume that, although there are 
different length scales present in reduced n-point correlation 
function, all of them have to be roughly of the same order 
to give a significant contribution. This is supported by the 
fact that - by its very construction - the reduced N-point 
correlation function vanishes when a single point or group 
of points from this set of N points are moved to large sep- 
aration (In that limit, the correlation function is just the 
product of lower order reduced correlation functions). For 
a geometrical picture, one can think of a polyhedron, in- 
scribed in sphere, with particles at each vertex having their 
velocities directed towards the centre of the sphere. (This 
configuration has the relative velocity of particles directed 
along their relative separation and hence can satisfy the sta- 
ble clustering ansatz.) For such a configuration, the scale in 
the correlation function will be the radius of the sphere cir- 
cumscribing the polyhedron. If the correlation functions are 
described by a single scale, then a natural generalisation of 
equation (4), will be 

In - {xf N ^)/x 3 ^ (5) 

The validity of such an ansatz is open to question and 
we do hope to check it directly in numerical simulations in 
a future work. In this paper we accept the above ansatz as 
a working hypothesis and use it to calculate the Sn param- 
eters in different regimes. Since these parameters have been 
studied numerically we can directly test predictions of this 
ansatz against existing results to obtain a feel for the valid- 
ity of the ansatz. It may be noted that, even though several 
models has been proposed to predict the values of Sn pa- 
rameters ( Hamilton 1988, Fry 1984, Schaeffer 1984, Balian 
& Schaeffer 1989, Bernardeau & Schaeffer 1992 ) they fail 
to reproduce correct values for lower order Sn parameters. 
We shall show that the non linear values of lower order Sn 
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parameters are predicted with fair degree of accuracy in our 
model. 



3 THE CORRELATION FUNCTIONS IN 
DIFFERENT REGIMES 

We shall now consider the implementation of the above 
ideas in three different regimes of gravitational clustering. 
We shall call the first one " perturbative regime " in which 
we expect perturbation theory is valid. The second regime 
( which we call " intermediate " regime ) is dominated by 
scale invariant radial infall of high peaks. Finally the third 
regime (" nonlinear ") is dominated by virialised blobs of 
matter. While we are mainly interested in the latter two 
regimes, we shall begin with some important observations 
regarding the perturbative regime. 

(a) Perturbative Regime 

We divide the density field into two parts at each point 
with one part coming from spherical collapse ( which we 
call the " monopole " part ) and the rest of the contribution 
comes from higher order spherical harmonics, characterising 
shear, tide and nonlinear coupling between them: 

5(x,a) = 8 sph (x,a) + e(x,a) (6) 

It should be noted that one can not assume S(x) = 5 sp h(x) 
for each point. While this may be obvious from symmetry 
considerations, a more formal argument can be given along 
the following lines: Let us assume for a moment that we can 
set e = 0. Since the growth of density contrast in spherical 
collapse model is well known (Peebles 1980), we can expand 
5sph{x,a) in taylor series to get 

oc oc 

6(x,a) = S sph (x,a) = ]T 5™{x,a) = £ ^<5 (1)Ar (7) 

JV=1 N=l 

where ^ 2 = 34/21, fi 3 = 682/189, m = 446, 440/43, 659.... It 
is clear that, in spherical collapse, S^ N ' oc S N . On the other 
hand, since (8) = we demand (8 (N ') = at every order 
of perturbation. This implies, in case of spherical collapse 
model, that (S N ) = for all N i.e., moments of S N vanish at 
all orders; hence we must conclude that 5 vanishes at each 
point identically. Clearly, we cannot set e = in (6) . 

For a generic Gaussian field, we have to work with a 
5 which has two parts, one coming from pure spherical col- 
lapse the other part, e is related to deviation from spherical 
collapse dynamics. The taylor series will then be 

oc oo oo 

5(x,a) = £ 5^(x,a) = £ ^ N + £ e N a N (8) 

JV=1 JV=1 ' N=2 

where we have expanded e in a perturbative series with ejv 
being of order 5 . ( Note that there is no contribution to 
e in linear order ). Although both the terms will be impor- 
tant for a generic random field, it has already been shown ( 
Bernardeau 1994a) that e becomes less dominating for rare 
events i.e. for large values of v —(8/a). In the perturbative 
regime, where a is small, any deviation from homogeneity is 
a rare event and hence one can assume that, statistically, e 
will be close to zero at most of the points. One can explic- 
itly demonstrate this claim by calculating the parameters 
Sn = (S N )c/(8 2 )i N ~ 1 ^ in the limit eat — > and showing 



that it will reproduce the well known results of Sn derived 
earlier by summing up all tree level diagrams in the limit 
o~ — ► (Bernardeau 1992 ). Consider, for example, the case 
of ft; we have 

S S = (8 3 )/(S 2 ) 2 = 3<5 (1) V 2) }/<5 (1)2 } 2 (9) 

= 3((5 (1)2 )(l/2/i 2 5 (1)2 + e2 ) 

+ M 1)2 )+M2(5 (1)2 ) 2 )/(5 (1) V) (10) 

The first term vanishes because (S^) = and the second 
term gives us vanishing contribution in the limit e 2 — > and 
we get the well known result S3 = 3/i 2 . A similar calcula- 
tion for higher order moments reproduce tree level results 
of perturbation theory, Si = 4^13 + 12/u| etc. These are the 
exact values of Sn parameters in the limit a 2 — > (i.e. at 
the tree level of perturbative calculation neglecting all loop 
corrections ) obtained previously by rigorous analytical cal- 
culations ( Bernardeau 1992, 1994a, 1994b, 1994c, 1995 ). 
Our analysis reconfirms that any deviation from spherical 
dynamics does not alter the values of Sn parameters at tree 
level in which only monopole part of the dynamics is rele- 
vant. The higher order harmonics ( shear, stress and their 
couplings ) start contributing only from loop level. ( This is 
also true for approximation schemes like Zeldovich approxi- 
mation etc; see Munshi et al. 1994 ) 

(b) Intermediate regime 

In the intermediate regime, we concentrate on the col- 
lapse of regions around peaks in which the density contrast 
scales as the correlation function. [We shall work with a 
fi = 1 universe]. Consider a spherical region of initial ra- 
dius Xi and overdensity S = v 2 ^L{xi) = a 2 v 2 x i '™ p+3 ' where 
n p is the index of the initial power spectrum and 00 is 
a constant. This region will expand to a maximum radius 
xta = (xi/S) oc v~ 2 x ( i np+4 ^ and then collapse back to a final 
radius Xf oc x ta - In the scale invariant radial collapse, the 
resulting profile will scale with x ta ( Fillmore & Goldreich 
1984, Bertschinger 1985, Hoffman & Shaham 1985). Taking 
x — \x ta and using equation (4), it is easy to see that 

£„(*) ~ (-gL) V )a; -3(" P +3)/(n P +4) (11) 

where we have introduced the notation z = 6/ (n p + 4). 

Evaluating the average < .... > using the Gaussian dis- 
tribution we find that the final result can be written as 
£2 (a:) = A^ lin (l) where I w x^ /3 and 

A ^=F(^ 2< *"" /2r * + H 6/ * (12) 

The above result is the generalisation [in the intermediate 
regime] of the analysis presented in Padmanabhan ( 1996 ) 
taking into account the averaging over different va peaks. It 
shows that the averaging introduces a spectrum dependent 
scaling. 

Let us now consider the higher moments. Using our 
ansatz for higher order moments (equation (5)) we can now 
compute the result for fjv to be 

I n{x) = (^2_) Z ^- 1 > { ^(iV-l) );c -3(iV-l)(n p+ 3)/(n p+ 4) ^ 
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The scaling we get for higher order moments is clearly hi- 
erarchical in nature. Using the definition of Sn parameters 
we find that in this intermediate regime 

sT = £n/& n - 1) = ( V (N - 1)Z )/( V Z ) (N - 1) 



or, equivalently, 



ST = (4tt) 



(JV-2J/2 



' z(N-l)+l ' 



(14) 



(15) 



we can a 

= S^T^ (^) = SvA^e^ (0 ( 16 ) 



Using the above results, we can also directly relate the £jv(a:) 
with £2(0 and otain 



fc,l Nonlinear Regime 

In this case, we take the initial density contrast to scale 
as the variance of the Gaussian random field, so that S oc 
x i ( n p+ 3 )/ 2 \y e assume as before that a patch with initial 
radius Xi will attain a maximum radius x ta = x/5 which 
will the collapse to form structure of size x = Xx ta - Then, a 
corresponding calculation gives 



£ a(a! ) = (^)» (l/ » )a .-3(»p+3)/(np+5) 



(17) 



where we have introduced the notation y — 6/(n p +5). After 
averaging over the initial Gaussian distribution, this result 
becomes £2 (a:) = B^2,iin(l) 3 ^ 2 where 



(v») 



y\3/y 



A 3 



6/y 



(18) 



This generalises the corresponding result of Padmanabhan 
(1996) to the nonlinear regime by taking into account the ini- 
tial Gaussian fluctuations. The averaging introduces a spec- 
trum dependent prefactor. 

For higher order moments analysis can be done in a 
equivalent way and the result is 



ffO\y(N-l) „(JV-1), -3(JV-l)(n p +4)/(n„+5) 



F I \ fU0\V("-i) 1 y(N-l) s 



(19) 



where y = 6/(n + 5). So the scaling we get for higher order 
moments is again hierarchical in nature and the Sn param- 
eters can be evaluated exactly in a same manner as before, 
giving: 



s N ° n = weT" 15 = (^ N - i)y )/(^) {N - l '> 

or, equivalently 



ST = (4tt) 



r ( y(N-l)+l \ 
(N-2)/2 L \ 2 ) 



(N-l) 



(20) 



(21) 



(22) 



This result can also be expressed as 

The averaging process < ... > in both quasi-linear and 
nonlinear regimes can be made more sophisticated by intro- 
ducing an additional weight factor which is proportional to 
some power of Lagrangian volume of the patch from which 
the object is collapsing i.e. X™ (see Padmanabhan et al. 
1996). In that case the results generalise to 



£n(x) = {x 



3(JV-l)+ro 



)/2{xT)x 



m\ 3(JV-1) 



(23) 



which can be simplified to 

S N = {2 (^}) N - 2 (^ N - 1 ^+ m ^/(p^+ m ^ ( - N -^ (24) 

where /3 = x/3 in quasi- linear regime and y/3 in nonlinear 
regime. Finally the Sn parameters are recovered after doing 
the averaging as before 



v-| l + „, <+i , 



(25) 



The simplest choice is m = which we shall use in this 
paper. It may be noted that the model used by Jain et al. 
(1995) corresponds to m — 3; the expressions given above 
can be used to read off the Sat parameters in any other 
scheme. (As we shall see in the next section, m = seem to 
give fairly good fit to the numerical simulation results ) . We 
may also note that: 

(i) The expressions derived for Sn parameters are valid 
for N > 2 (Si = 1 by definition ). 

(ii) Direct comparison with results of intermediate and 
nonlinear regime shows STi n p) = S N on (n p + 1). Also note 
that the value of Sn is independent of A. 

(iii) Temporal dependance of £jv in both quasi-linear 
and nonlinear regime can be derived from the fact that any 
statistic of scale invariant system can be expressed as a func- 
tion of x/xni, where x n i is the scale defined through the re- 
lation a(x„i) = 1. Since x n i oc a 2 /( n + 3 ) all correlations will 
be function of single variable q = xa~ < - n+3 ^ 2 . 

(iv) It is clear that except for calculating the averages of 
powers of v ( which is assumed to be distributed normally 
) nowhere have we actually used the fact that the initial 
density distribution was Gaussian, which clearly show that 
our method of analysis can be generalised in a straight for- 
ward manner to calculate Sn parameters for initially non 
Gaussian distributions. 

(v) For studying gravitational clustering in dimensions 
other than 3 the same method of analysis can be used with 



the scaling ^n(x,o) oc £, 
and £n(x, a) oc £ 



D(JV-l) 



2 ./in 



D(N-l)/2 



(I, a) in intermediate regime 
(I, a) in highly nonlinear regime. 



(vi) Given the Sn parameters, one can compute the void 
probability distribution function and related quantities. This 
calculation is indicated in the Appendix. 

(d) Transition between the regimes 

Having determined the behaviour of correlation func- 
tions in the three different regimes, one can enquire where 
the transition between the regimes occur. Since there ex- 
ists three distinct phases in gravitational clustering we have 
two transition points: (1) Transition from the perturbative 
regime to intermediate regime and (2) Transition from in- 
termediate regime to nonlinear regime. Let the first tran- 
sition occur when £,2,Un{l) = and the second when 
^2,(in(0 = Tc 2 ■ Finding the value of ^2,i»n(2) for which quasi- 
linear and intermediate £2(2) matches we get Tj^ = 1/A 1 ^ 2 . 
Similarly, equating the expressions for the intermediate and 
nonlinear regimes gives = (B/A) 2 ^ z . 

We have used the two point correlation function to de- 
fine the transitions since they are most directly related to 
the density inhomogeneity. It is, of course, possible to repeat 
the same exercise using higher order correlation functions. 
Our results for the higher order correlation functions can be 
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Figure 1. Transition points T C1 and T C2 are plotted as a func- 
tion of spectral index n p . The short-dash curve corresponds to 
transition points from perturbative regime to intermediate regime 
while the long-dash correspond to transition from intermediate 
to nonlinear regime, predicted from our model. The curves are 
normalised to match Jain et al.'s simulation result for n = —1. 
Circles corresponds to values obtained in the simulation by Jain 
et al. (1995). Transition points were obtained by finding the inter- 
section of straight lines which we fit to represent various regimes. 



the two rgimes together by using a simple approximation. 
We shall discuss this approach in this section. The results 
of last section can be obtained as a special case of this ap- 
proach. To do this we begin with equation (H) written as 



3h 



(29) 



dP _ dP 

OA ~ dX 

where we have introduced the following new variables 

D = ln[(l + f2(a;,a)], A = In a, X = \nx (30) 

and written the pair velocity as v = —hax. Simulations in- 
dicate that h m 2 in intermediate regime and and h ~ 1 in 
nonlinear regime which we have used for results obtained in 
earlier sections. Here we will try to get an unified picture 
covering both the regimes. We shall now assume that we 
can treat h as approximately constant while integrating this 
equation. In that case, the general solution is 



1 + 6 = a 3h F(a h x) 



(31) 



where F is an arbitrary function to be determined by initial 
condition. If the linear £l is a power law, we know tht the 
true £2(0, x) can only depend on the variable q = xa~ 2 ^ n+3 ^ 
which is possible only if F is a powerlaw. So we must have 



£2(3;, a) oc a (ax ) 



(32) 



The index 7 can be determined by matching the above ex- 
pression with linear two point correlation function at a scale 
x c = a 2 /(™+ 3 ) which is going nonlinear. We then get 

3h(n + 3) 



7= M- + 3)+2 (33) 
Now it is possible to write the two point correlation function 



summarised as 
(x) 



f- n on 
SAT 



Qtree — 1 /j\ 

1 c non i3(JV-l)/2 , n 



a N— 1 r,int j- 



B 



(26) 



Using these, it is easy to see that the transition points de- 
fined through N-point correlation function will give us 



T. 



(JV) 



/ rttree 

— WN 



is. 



int\l/2(N-l)rp(2) 



N 



and 

T (N) 



/ Qnon 



I Q int\2/3(N-l)rp{2) 
I^N ) 1 c 2 



(27) 



(28) 



Since S% ee < S%* and S'^ > S% m , it is clear that tran- 
sition for higher order moments will occur for smaller and 
smaller values of £2(1)- 

It should be noted that all though Sn parameters are 
insensitive to the modelling parameters like A the transition 
points are sensitive to the choice of these variables. We have 
taken A = 1/2 which is close to value taken by Jain et al. 
(1995) for their fitting function. 



4 UNIFIED ANALYSIS FOR DIFFERENT 
REGIMES 

In the last section, we discussed the intermediate and non- 
linear regimes separately. It is, however, possible to discuss 



£ 2 ( X ,a) (X a eh /l 2 + h ( n + 3 )] x - 3h (n + 3)/[2 + h(n + 3)] 



(34) 



Earlier results of intermediate and highly non-linear regime 
can be recovered by taking h — 2 and h = 1 respectively. 

It is actually possible to relax the power law require- 
ment and still obtain a general result. From the character- 
istics of equation ( |29| ) we can show that £2(2,0) can be ex- 
pressed as a function of a) where I 3 = x 3 [l + £2(2;, a )]- 
That is 

Za{x,a) = U[Zz,(l,a)] (35) 

where U is some function. Combining with ( |3l| ) we have 
a 3h F(a h x) = U[£ L (l,a)] or equivalently a 3h F(r) = U[a 2 Q{r)} 
where we have used the fact that we can write I 3 = x 3 fe = 
r 3 F(r) with r — a h x. But the above expression can be valid 
for arbitrary a at fixed r, only if 

U(z) oc z 3h/2 (36) 

So, in terms of correlations functions, we must have 

f a (z,a)«[&(!,o)] ah / a (37) 

This generalises the relations £2 oc 1^, £2 oc t^ 2 we used for 
quasilinear and nonlinear regimes in the last section. 

To perform the averages over regions with different peak 
heights, we only have to do the rescaling fx, — ■> v 2 ^l and 
note that £2 oc (xt/x) 3 oc £^ 2 (x») implies x oc Xi^/ 2 (xi). 
Simple algebra then gives 

6 ( X)a ) oc <fi> oc (l ,6V(2+M™+3)) )a .-3Mn+3)/(2+M™+3)) (gg) 
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This generalises the relation ( [tl"| ) and (^) of the last section. 
Assuming that v is a Gaussian variable and performing the 
average, we get upto a normalization, 



6(x,a)oc2^ 2 r(^) 



r . O-j-1 , _3fc(„ + 8)/(2+M»+3)) 



(39) 



where q = 6/i/[2 + fa(n + 3)]. Normalizing the expression 
properly, we can write the final result as 



t 2 {a,x) = A{h,n)[&{l,a)f h/2 
where 

'2x¥rr(^i) 



A(h, n) 



(I)' 



2^ 



(40) 



(41) 



To obtain the earlier results for intermediate regime we have 
to set h = 2 which gives q = 6/(n + 4) and 



.4 













- 2v^T J 



(42) 



similarly for the nonlinear regime we have h = 1, a = 6/(71+ 
5) and 



B = 



2^3/2 

A 



v 2 , 



20F 



(43) 



These results match with earlier expressions. 

Using our ansatz in (JE]) it is possible to generalize the 
result for higher order correlation functions. We find 



which can be explicitely written as 



(44) 



I 2 

r(s±i) 



(45) 



This allows calculation of Sn((,2(x, a,)), given hfaix, a)). 



5 COMPARISON WITH SIMULATIONS 

The results obtained in the previous section can be com- 
pared with the simulations as regards two essential aspects. 

First of all, the results show that the spectral depen- 
dence of the scaling relations between the nonlinear and lin- 
ear correlation functions can arise due to averaging peaks 
of different heights. This, in turn, implies that the scales at 
which the transition from perturbative regime to intermedi- 
ate regime, or from intermediate regime to nonlinear regime 
takes place depends on the power spectrum index. By com- 
paring the predicted values for these transitions with the 
result of simulations, we can test the the validity of our 
averaging scheme and the basic model for the two point 
correlation function. Secondly, we can compare the values 
of £n obtained from the model with those of simulation in 
both intermediate and nonlinear regimes. Since £2 is related 
to ^2L in a nonlocal manner, and our ansatz relates £jv to 
£2 (in an indirect manner), we would expect some nonlocal 
relationship between £jv and £2- This will test the validity 
of our ansatz regarding higher order correlation functions. 
Note that these two comparisons test the two distinct gen- 
eralisations of the work in Padmanabhan (1996), introduced 
in this paper. 

The comparison of predicted values for the transition 



is shown in figure 1 along with results of numerical simula- 
tion given in Jain et al. (1995). We see that there is good 
agreement between theory and simulations suggesting that 
(i) the basic picture for the evolution of gravitational clus- 
tering, developed in Padmanabhan (1996) is correct and (ii) 
the spectrum depedence of the scaling relation can be under- 
stood by averaging over the initial fluctuations. The analysis 
also shows that — as n p changes from -2 to — the lower 
transition point varies between 0.25 and 1.0 while the upper 
one varies between 2.0 and 7.0 

We shall next turn to the comparison of £jv predicted 
by our model with the results of numerical simulations, in 
order to test the validity of our ansatz. In doing so, one 
should be aware of several effects which could "corrupt" the 
values of Sn parameters in the simulations, and take ade- 
quate precautions to correct for them. At small a, i.e. in the 
perturbative regime, the main contribution to error comes 
from cosmic variance i.e. due to presence of small number of 
large cells containing completely independent samples. This 
error can be reduced only by increasing the size of N-body 
computation box. On the other hand in the highly nonlinear 
regime one is restricted by resolution of the N-body simu- 
lation for probing very small scale. Poisson shot noise also 
starts playing increasingly dominant role as soon as the av- 
erage occupancy of cells becomes comparable to unity. The 
usual procedure used for computing Sn parameters is by 
taking moments of cell counts P(n) for different cell sizes. 
In general, cell counts show a power law behaviour in the 
highly non linear regime upto n — n c (— n^) followed by 
an exponential tail at n > n c (Balian & Schaeffer 89). In 
an ideal (infinite) catalogue this exponential tail will be ex- 
tended to very small values of P(n); but, in practice, there 
is a sharp cutoff around n — n max which is the most dense 
cell present in the N-body catalogue. Higher moments of 
P(n) - and hence higher Sn parameters - are more sensitive 
to large n tail of P(n). So it is extremely important to ex- 
tend the the exponential tail to infinity and then normalise 
the corrected P(n) again before calculating Sn parameter. 
This way of correcting for measured Sn has been extensively 
studied and expected to give correct - or at least, more re- 
liable - result (Colombi et al. 92, 94, 95, Lucchin et al. 94 ). 
Recently Colombi et al. ( 1995 ) has done a careful analy- 
sis of high resolution n-body data with large dynamic range 
correcting for all the errors mentioned above. Their study 
covers power law models n v = 1, 0, -1, -2. and they study 
first three non-trivial Sn parameters i.e. 5*3, Si, and S5. 

We have computed |jv for N = 2,3,4,5 for the power 
laws n = —2,-1,0, 1 using the published data in Colombi 
etal. (1995) and our ansatz. The results are shown in fig- 
ures 2 to 5. The theoretical prediction based on our ansatz 
is shown by solid, S-shaped line and the simulation data 
is shown by points. The straight lines have the asymptotic 
theoretical slopes. 

It is clear from the graphs that our basic ansatz is an 
attempt in the right direction. The overall agreement be- 
tween the theory and simulations is good especially when 
we consider the simplicity of our model. We will now com- 
ment on several details related to the comparison between 
theory and simulations. 

Plots of £n(x,o) vs £2 (x, a) shows that our basic claim 
regarding three phases in graviational clustering seems to be 
correct and also one do get a nonlocal scaling relation for 
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n = -2 




Figure 2. £is[(x,a) has been plotted against £2(^1) for to = — 2 spectra. Solid curve in panel S2 corresponds to fit by Jain et. al. ( 
1995) and dash curve for Hamilton et. al. ( 1991), in other panels solid curves correspond to prediction from our model. Straight lines 
in different panel correspons to slopes (TV — 1), 3(7V — 1) and 3(N — l)/2 for perturbative, intermediate and highly nonlinear regime 
respectively. Dots represent N-Body simulation data from Colombi et. al. ( 1996 ) 




Figure 3. Same as figure 3 for n = —1 spectrum 
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Figure 5. Same as figure 3 for n = 1 spectrum 



© 0000 RAS, MNRAS 000, 000,000 



hivolution oj correlation junctions 9 



N-point correlation function similar to scaling for two point 
correlation function as suggested from the characteristics. To 
some extent the scaling in higher order correlation function 
reflect underlying scaling in two point correlation function; 
but it can also be argued that if £ N (x,a) can be expressed 
as a smoothly varying function of £2(0:, a), £n(x,o) = 
Tn(&(x, a)) then we can write £n(x, a) = TN(F n (£2(l, a))) = 
G u ,n(£~2(1, a)) where G nj jv = Tjv * F n and thus define a scal- 
ing relation between £n(x,o) and £2(1,0). 

The relation between £2(0,2:) and £2(1,(1) shows good 
agreement with the analytical fit suggested by Jain et. al. 
(1995) The scatter in the data increases with n, which can 
be understood in following manner. We use x and £2(1, a) 
to recover the lagrangian radius I = x(l + £~2(x, a)) 1 which 
was then used to get the linearly extraploated £2(1,0) = 
a 2 a 2 l~ l ' n+3 \ Error in estimation of x or £~2(x,a) from the 
published data of Colombi et al. (1995) gets reflected finaly 
in error of £2(1,0,). We can relate fractional error A^ 2( - ; a ) 
with fractional error Aj by |Af 2 (; a )| = (ra + 3)|A;| Which 
shows that for same A;, A^^ j increase with n. This ex- 
plains (partly) why we get more scatter for n = 1 spectra 
compared to n = — 2 spectra. 
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APPENDIX - A 

Given the Sn parameters, one can compute the void prob- 
ability function which is of some theoretical and practical 
importance. In scaling models (which assume Sn parame- 
ters are constant over some length scale) the void probability 
function can be written as ( White 1979 ) 



P = exp(-0(n c )/£ 2 ) 



(46) 



where <j> is generating function for Sn parameters and is 
defined as 



(47) 



where n c = n£2 is a scaling variable. Substituting our ex- 
pression for Sn in the above equation we get 



(j)(n c ) =n c - y^(47r; 



( P -2)/ 2 r(((p-l)s + l)/2) 

r(( s + i)/ 2 ))p-V! 



(-nc) p (48) 
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where s = z for intermediate regime and s = y for nonlinear 
regime. Using definition of gamma function it is easy to write 
down the sum as 



K) = ,-2<i±M x 



00 ,. c 

£ / 

v =2 J 



47T 

dtt (p ~ 1)s/2 exp(-t) 

Vtp\ 



n(s + i)mj 



(49) 



Interchanging the sum and the integral we get 

^ ^ _„ r[( 8 + i)/2] v 



f 

Jo 



dtt 



47T 

-(s+l)/2 



e-*(exp(-/t s/2 ) + /t 



s/2 



where 
/ = 



47T n c 



r[(* + i)/2] 



1)(50) 



(51) 



and taking the limit n c — > 00 shows that 4>(n c ) w n c /2 
for large n c . In the scaling model proposed by Balian and 
Schaeffer, <j> is expected to scale as n]r^ asymptotically; 
hence our ansatz leads to the uj = model. (Note that other 
extreme case is ui — 1 which is the negative binomial model 
proposed earlier). The counts-in-cell P(n) is proportional to 
n" -2 showing that, in our model, P(n) oc rT 2 . 
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